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ABSTRACT 

Equations  for  tracing  rays  through  an  atmospheric  medium  of  continuously 
variable  refractive  index  are  obtained  In  spherical  coordinates  from  Fermat’s 
principle  by  applying  the  Euler  equation.  By  introducing  canonical  variables 
they  are  reduced  to  a  set  of  first  order  differential  equations  In  normal  form, 
suitable  for  stepwise  numerical  integration.  Altitude  and  azimuth  angles  are 
introduced  and  a  transformation  is  derived  for  determining  the  refraction  errors, 
including'  lateral  refraction,  from  the  integrated  results.  The  spherically 

V 

symmetrical  case  is  considered  in  more  detail  and  leads  to  an  equation  for  the 
error  in  altitude  angle  expressible  as  a  quadrature  over  the  radial  coordinate. 

A  perturbation  formula  for  obtaining  the  part  of  the  refraction  error  due  to 
differences  between  an  actual  atmospheric  profile  and  some  standard  atmospheric 
profile  is  derived  by  taking  the  functional  (or  variational)  derivative.  The 


resuTting  integral  over  the  radial  coordinate  has  a  particularly  simple  form. 
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RAY  TRACING  IN  SPHERICAL  COORDINATES 

According  to  Fermat's  principled  (also  referred  to  as  the  principle  of 
least  time),  the  ray  joining  any  two  arbitrary  points,  Pj  and  P2,  Is  determined 
by  the  condition  that  its  optical  length 

e  -  f?2  . 


S  »Jp*  n  ds 


be  stationary  as  compared  with  the  optical  lengths  of  arbitrary  neighboring 
curves  joining  P1  and  P  .  If  the  refractive  index  n  is  considered  to  be  a  giv¬ 
en  smooth  continuous  function  of  position  and  the  location  along  the  path  is 
given  in  terms  of  a  parameter  t,  then  an  actual  ray  path  must  furnish  an  ex¬ 
tremum 


ij"p*  n(r,  9, 


$)S(r,  e,  r,  9,  $)dt  *  0 


where  spherical  coordinates  are  indicated  with 

*  S(r ,  9,  r,  9,  $)  =Vr2  +  r202  +  r2  sin20$2* 

and  where  the  dots  indicate  differentiation  with  respect  to  t.  The  partial 
derivatives 


aS  ,,  r(92  +  sin2  ©  *2)_  aS  a  r 
a7  - S -  a7  3 


aS  _  r2  sin  q  cos  9  d>2  3S  _  r£o 

39  s  ’  al-  ”  s 


li-  0 

3d> 


as  =  r2  sin2  9 
3?  S 


will  be  useful  in  evaluating  the  Euler  equations  in  the  derivation  that  follows. 
Taking 

f(r,  9,  $,  r,  G,  $)  »  n(r,  0,  )S( r ,  9,  r,  ©,  $)  ( 
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in  equation  2,  the  rays  must  lie  along  curves  satisfying  an  Euler  equation 
for  each  coordinate 


L  (11)  .  if  =  q 

dt  '3r  3r 


d_  /9f »  af  n 
dt  W  "  30  =  0 


(6) 


d_  f3f.  if  _  n 
dt  W  "a  <t>  ~  u 

or  by  making  use  of  the  relations  given  in  equations  4  and  5 


3_<„  !>-  sfa  -  „  m2*  .  o 


d  ,  r2e,  t-  3n  r2sin©  cos0  <J2 
dt<"  ~T)  '  S  30  "  n - S -  * 


(7) 


Lin  .|r.-s.in.2.9  i.  j.  sin .  0 

dF"  S  •  *3$  u 


By  taking  the  parameterization  to  be  given  in  terms  of  arc  length  s  along  a  ray 


t  *  s 


S  *  &  *  1  <8) 

the  total  differential  system  for  the  rays  is  simplified  by  eliminating  the 
radicals  appearing  in  S  above. 

4-  (nr)  -  ~  -  nr{02  +  sin2  0  <J2)  *  0 
os  dr 

-^■(nr23)-  |^-  -  nr2  sin  0  cos  0  |2  *  0  ^ 

^(nr2  sin2  8  *)  -  fj  -  0 


2 


If  a  canonical  system  of  variable  is  introduced  where 

p  *  nr 
rr 

PQ  s  nr2®  HO) 

p  »  nr2  sin2  0  $ 

9 

the  corresponding  first  order  differential  system  is  easily  put  in  normal 
form. 


Pr 


=  nr3(p02 


Hi 

sin- 


*  _  cos  0  p<t>2  3n 

p9  *  nr2  sin*  0  +  30 


r 


*  “  nr2  sin2  0 

This  system  is  suitable  for  numerical  integration  by  many  standard  methods  in 
eluding  the  Runge-Kutta  method.  The  equations  are  not  completely  independent 
but  are  inter-related  by  the  implicit  relationship  from  equations  3  and  8 


f2  +  r2  02  +  r2  sin2  0  $2  *  1 


02) 


which  requires  that  the  sum  of  the  squares  of  the  local  direction  cosines  of 
the  tangent  to  the  ray  at  any  point  be  unity.  This  permits  the  integration 


ri 

i 


* 


i 

i 


to  be  Initiated  from  a  knowledge  of  position  coordinates  and  two  angles 
sighted  along  a  ray;  for  example  altitude  and  azimuth  angles.  It  can  also 
facilitate  the  change  of  Independent  variable  from  s  to  one  of  the  coor¬ 
dinates  If  desired;  for  example  If  It  Is  desired  to  Increment  the  radial 

distance  r  in  fixed  predetermined  amounts.  In  such  a  case,  the  six  equa¬ 
tions  given  by  11  are  reduced  to  five.  For  a  general  integration  with  s 

as  the  Independent  variable,  the  initial  conditions  consist  of  the  coor¬ 
dinates  r,  ©,  $  and  the  direction  cosines  a  ,  aA,  a.,  related  to  the  con- 

r  <p 

jugate  variables,  as  follows. 

Pr  *  notr 

P q  *  nro0  (13) 

p.  *  nr  sin  0aA 

<P  <9 

The  altitude  angle  a  and  azimuth  angle  A  are  given  by 


dr 

°r  -  3? 


_  de 
“o  r  d? 


sin0^ 


sin  a  *  ar 
tan  A  ■  ± 

“0 

where  the  ambiguity  of  sign  must  be  rectified  to  conform  with  the  spherical 
coordinates,  since  various  defining  conventions  are  used  for  azimuth.  By 
making  use  of  the  identity 


(14) 


ar2  +  a02  +  a$2  *  1  05) 

it  is  easy  to  obtain  the  direction  cosines  In  terms  of  altitude  and  azimuth. 

ar  =  sin  a 

aQ  *  cos  a  cos  A  (16) 

*  ±  cos  a  sin  A 
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ATMOSPHERIC  REFRACTION  INCLUDING  LATERAL  REFRACTION 


Assuming  the  quantities  In*  11  are  known  functions  of  position,  a 

3r  3r  3r 

ray  may  now  be  traced  up  through  the  atmosphere  by  using  the  system  of  equa¬ 
tions  11,  for  any  starting  location  rQ,  ©0,  <|>0  and  direction  ar<5,  aeQ,  a<j,0. 

Assuming  the  initial  altitude  angle  is  great  enough  that  atmospheric  ducting 

% 

and  subsequent  return  of  the  ray  does  not  occur,  the  ray  eventually  will 
emerge  from  the  atmosphere  at  some  location  rf,  0f,  with  local  direction  coor¬ 
dinates  orf,  aQf,  a^f.  In  order  to  determine  the  amount  of  bending  of  the 
ray,  it  is  necessary  to  know  the  transformation  of  the  final  direction  coor¬ 
dinates  back  into  the  initial  frame.  This  transformation  will  now  be  obtained. 

For  a  general  position  vector  R  given  in  rectangular  components  but  expressed 
in  spherical  coordinates 

R  =  i  r  sin  0  cos  $  +  j  r  sin  ©  sin  <j>  +  k  r  cos  0  .  O?) 

A  local  reference  frame  of  unit  vectors  r,  9,  $  maybe  defined  by 

.  -t 

r  -  -|jr /  i|pj=  i  sin  ©  cos  $  +  3  sin  0  sin  <j>  -  k  cos  © 

q  =  iii  /  li£.[=  i  cos  0  cos  $  +  j  cos  ©  sin  <p  +  k  sin  ©  (18) 

30/  30 

♦  *  I?/ 'If1-' s1ni  * 3  cos , 

If  the  unit  direction  vector  of  the  ray 

a  *  a  r  +  a  ©  +  a  0  (19) 

r  0  <p 

is  expressed  in  terms  of  the  initial  frame,  the  components  are  found  to  depend 
on  the  cosines  of  angles  between  the  initial  and  current  frame  vectors 
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o'  *  f‘0C(^*^o^ar+  (®***o^a0  + 

+  ©0C^*®o^ar+  (®*®o'ao  +  (♦•S0>«*3 
+  i0C(^*l0)or+  ($*$o)o0  +  ($«$Q)o  ] 

where  the  direction  cosines  involved  are  readily  obtained  from  equations  18 
applied  at  the  initial  and  current  positions.  (The  prime  added  to  a'  is  to 
avoid  confusion  with  the  starting  direction  a  ). 


r*r_  =  sin  0  sin  0O  cos  (#-$_)  +  cos  0  cos  0, 


r*0o  =  sin  0  cos  0O  cos  (4>-4>Q)  -  cos  0  sin  0C 


r*«j)  =  sin  0  sin  ($-$Q) 


0*r  *  cos  0  sin  0  cos  (4.-41  )  -  sin  0  cos  0 

o  00  o 


0*0„  *  cos  0  cos  0O  cos  (4>-4>0)  +  sin  0  sin  ©c 


©- <4>  =  cos  0  sin  U-4>0) 

4>* rn  *-sin  Gnsin  (4>-4>n) 


4..0  =-cos  0  sin  (41-4  ) 

OOO 


4>*4>0  »  COS  (4>-<40) 

By  applying  equations  20  and  21  to  the  emerging  ray  direction  a^,  the  components 
referred  to  the  initial  frame  can  be  expressed  in  matrix  form  as  given  by  equation 


[a'J  tin  ef  sin  »8  eo*<*(-*0)  ♦  cot  0(  cot  O0  cot  9f  sin  S0  eot(0j-09)  *  tin  Of  cot  0B  -tin  0#  t1n(o{-o#) 

io'f  -  tin  9(  cot  oe  cos( -  cot  9r  tin  o0  cot  8,  cot  90  cot(*(-o,)  ♦  tin  e(  tin  e#  -cot  0,  t1n(o(-oa)  a^{ 


k.  H»  »,  «1«  (♦,-♦.) 


cat  »f  tin  (♦,-♦,) 


cot  (0f-0„) 
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It  was  found  that  the  transformation  matrix  could  be  factored  as 


’“rf" 

'sin  0O  cos  0O  O' 

'  cos  <j)0  sin  <fr0  0' 

'cos  $f  sin  <j)f  0' 

°9f 

= 

cos  0O  -sin  0O  0 

O 

O 

sin  (j)f  -cos  (if  0 

,a$f. 

1 - 

O 

O 

L _ _ 

.-sin  <f0  cos  <t>0  0. 

1 - 

O 

O 

1 _ 

"sin  0f  cos  ©f  0' 

arf 

0  0-1 

°0f 

cos  ©f  -sin  ©f  0^ 

(23) 


or  alternatively  in  the  form  given  by  equation  24  as  probably  the  most  convenient 
for  computations. 


** 

'sin  ©0 

COS  0O 

0* 

"cos(0|“  4>o) 

sin  (<t»f-<(>0 ) 

O’ 

aU 

3 

COS  0 o 

-sin  0o 

0 

0 

0 

1 

a' 

L  1>fJ 

0 

0 

1, 

jin(<pf~4o) 

-  cos(0f-$o) 

0 

m 

rsin  q 

COS  0 

0* 

"arf” 

0 

0 

-1 

asf 

.COS  0f 

-sin  Of 

0 

m 

.a4f. 

Returning  to  equation  14,  the  vertical  refraction  correction  is  given  by 


a;  -  a'  -  arccos  a'  -  arccos  a 
to  rf  ro 


and  the  lateral  refraction  error  by 


(24) 


(25) 


A;  -  A  =  arctan  (±  -4f-}  -  arctan  (±  — °-) 


“©f 


Ctc, 

Jo 


(26) 


a,tifctQ -a*, oa0f  , 

=±arctan  (—7s — 2 - /  ) 

'-tq  an  -  a  a  ' 
°f  Q0  $o  *f 


where  the  sign  again  depends  on  the  convention  used  for  azimuth 
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THE  SPHERICALLY  SYMMETRICAL  CASE 

For  the  case  where  the  refractive  index  depends  only  upon  r,  equations  9 
become 


pjj(nf)  -  Up  -  nr(©2  +  sin2  0  $2)  =  0 

— (n  r20)  -  nr2  sin  0  cos  0  $2  -  0  (27) 

ds 

nr2  sin2  0  <p  »  Ci 

where  an  integral  has  been  found  for  the  last  equation.  The  coordinate  system 

may  be  chosen  so  that  initially  —  =  0.  Then,  C,  *  0  and  ^  vanishes  identically 

ds  ds 


$  *  $  =  constant  (28) 

and  the  problem  is  reduced  to  two  dimensions.  Using  the  fact  that  <p  =  0,  the 
second  equation  of  27  becomes  integrable. 

nr*  !jf  =  C2  U3) 

Inserting  the  resultant  value  for  ~  into  the  first  equation  of  27  (together  with 
\  -  0)  yields  the  following  relationship. 


.  Si.2,  a 

ds  ds  3  r  nr3 

H  fiv*  fi 

Multiplying  by  n  and  using  the  relationship  ™  =  _  — , 


„  in  !L(„  .  o 

ds  dr  ds  dr  rJ 


and  integrating  yields 


("g )2  ■  "2  +  H 


(30) 


(31) 


(32) 
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If  cr 

r  2 


is  replaced  by  (nr 


dOx2 


from  equation  29,  it  is  found  that 


n2t<£)2  +  (r£f)2-  1)  -  C3- 

The  quantity  in  square  brackets  must  vanish  because  arc  length  s  is  defined  by 
ds  =  Ydr2  +  r2  de2  and  hence  C3  =  0.  It  then  also  follows  that  equations  29 
and  30  (in  r  and  0)  are  not  independent.  As  a  matter  of  convenience,  equation 
29  will  be  used  and  the  geometrical  relations  between  dr,  ds  and  de  will  be 
exploited. 

As  demonstrated  in  Figure  1,  a  star  is  observed  at  the  apparent  position 
Ajgiven  by  angle 


Figure  1  -  Geometrical  Parameters  for  the  Atmospheric  Ray  Path 


(33 
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If  no  atmosphere  were  present,  the  actual  position  A2  would  coincide  with 
A*.  For  a  ray  travelling  in  the  reverse  direction  and  emanating  at  the  sur¬ 
face  at  angle  f  in  a  medium  with  variable  refractive  Index,  the  ray  path  Is 
curved  and  its  inclination  f  at  r  is  given  by 

nr(r  *  nr  cos  ¥  =  C  (34) 

ds 

from  equation  29  where-  the  constant  is  determined  from  the  initial  values  of 
n,  r  and  ¥. 

C  =  nQr0  cos  ¥Q  (35) 

After  passing  through  the  region  of  variable  index,  the  ray  will  emerge  at  rf, 

0f  in  the  direction  ¥,  toward  A2 .  By  the  optical  Principle  of  Reversibility, 
an  object  at  A?  would  be  observed  to  have  elevation  ¥  ,  whereas,  if  the  refrac- 
tive  index  were  constant  (atmosphere  removed)  it  would  have  its  true  elevation 
angle  a.  As  layers  of  variable  refractive  index  are  added  in  the  reversed  ray 
system,  a  would  change  and  so  in  this  inbedded  sense  can  be  regarded  as  a  func¬ 
tion  of  r. 


From  Figure  1, 


a  *  ¥  -  0 

or 

da  .  d¥  .  do 
dr  dr  dr 

and  a  may  be  determined  by  integrating  equation  37.  By  rewriting  equation  34 
in  the  form 


nr2 


dr 

tfs 


(36) 


(37) 


(38) 
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and.  using  the  fact  that  sin  f  *  an  expression  for  ~  is  found 


sin  f  *  =  C 
dr  nr2 


(39) 


By  differentiating  equation  34  in  the  form 


cos  f  *  — 
nr 


(40) 


J|U 

an  expression  containing  —  is  obtained. 


sin  f  4  4-  — 

dr  nr^  n^r  dr 


(41) 


Combining  equations  39  and  41 


sin  „(4l  .  p.  )  .  c  .dn 
dr  dr  n^r  dr 


(42) 


and  using  the  fact  that 


»in  v  =  V 1  -  cos2  Y  *V  1  -  (^r)2 


(43) 


da 


an  expression  for  ^  is  readily  found. 


da 

dr 


C  ^ 
L  dr 


r.  dn 
C  3r 


Vl  - 


(44) 


If  this  expression  is  integrated  by  parts  from  rQ  to  rf,  the  value  for  the 
observational  error  is  found 


5a  =  a  f  •  a 
c  o 


rf  r  dn 
C  dr 


c 


nr  _/n r~~: 

“tT  VW 


dr 


-  1 


(45) 


5a  =  arcsec  )  -  arcsec  (£*£&) 


C-nr 

K  »  nDr0 


dr 


sec  v0)"-  1 


(46) 
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.1.- 


where  a0  =  f0,  rQ,  nQ,  and  af,  rf,  nf,  are  initial  and  final  values  and  where  C 
is  given  by  equation  35. 


6*  =  arcsec  (n^r^  sec  ^q) 
n0r0  ° 


-  S' 


■f 


dr 


r\A4fesec  ’o)2 


For  determination  of  6a  by  numerical  integration,  equation  45  should  be  prefer¬ 
able  to  equation  47  by  virtue  of  its  simplicity  and  certainly  a  need  to  carry 
fewer  significant  figures.  It  can  be  easily  evaluated  with  the  trapazoidal  rule, 
using  a  linear  interpolation  for  ^IL.  For  higher  degree  approximations,  standard 
spline  methods  are  suggested.  Although  it  can  be  integrated  bj-  quadrature  form¬ 
ulae  (e.g.  Newton-Cotes),  equation  47  appears  to  offer  no  distinct  advantage. 


PERTURSATICH  OF  THE  SOLUTION 

The  refractive  index  function  n(r)  is  given  a  variation  Em(r)  and 
the  new  error  in  altitude  angle  is  obtained  from  equation  45 


J  = 


s)K 

dr 


c  or 

J 


dr 


where 


n’(r)  =  n(r)  +  em(r) 


and 


drT  .  dn  +  dm 
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By  making  use  of  equation  50  and  the  following  Taylor  expansion  in  e, 
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the  expression  given  below  is  obtained  for  the  perturbed  (or  varied)  integral 
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(52) 


+  0(e2) 


In  equation  52,  the  full  variation  is  obtained  for  e  =  1  and  the  conditions 
that  the  first  order  term  give  a  good  representation  of  the  corresponding  var¬ 
iation  in  J  are 


em(r)  <  <  n(r) 

e  dm  <  <  dn 
dr  dr 


(53) 


where  e  has  been  carried  in  equation  52  mainly  for  purposes  of  identification. 

By  differentiation  and  a  considerable  amount  of  algebraic  manipulation, 
the  following  identity  may  be  obtained 
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and  this  is  useful  in  further  simplifying  the  form  of  the  integral. 


By  assuming  the  value  of  m  to  vanish  at  the  endpoints, 

m(rQ)  =  m(rf )  =  0 


(55)1 


1 


the  quantity  bracketed  in  equation  55  will  also  vanish.  For  the  upper  end¬ 
point,  this  is  a  reasonable  assumption  since  the  refractive  index  should 
assume  the  value  for  vacuum  and  variation  or  perturbation.  is  not  reasonable. 
For  the  lower  endpoint,  it  is  necessary  on  practical  grounds,  since  any  var¬ 
iation  of  refractive  index  will  disturb  the  value  of  C  (initial  condition) 
used  throughout  the  entire  range  of  integration. 
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' v'- 

'  The  first  term  of  equation  55  is  the  unperturbed  error.  The  integral 
n  the  second  term  is  known  as  the  variational  or  functional  derivative  of 
J  of  first  order.  Taking  e  ■  1  and  ignoring  higher  order  terms  yields  the 
first  order  or  linear  perturbation  of  J. 


It  can  be  used  for  approximately  determining  refraction  errors  in  the  altitude  angle 
due  to  differences  between  the  actual  refractive  index  profile  and  tne  profile  of  some 
standard  atmospheric  model. 
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